%% Load the data - Select the File: 'BulkViscosity'
close all; clear all; clc;
clearvars;
total_directory=uigetdir();
total_directory = dir(total_directory);
dirFlags = [total_directory.isdir];
total_directory = total_directory(dirFlags);
total_directory = total_directory(~ismember({total_directory.name},{'.','..'}));
'Done Load'
% Set variables for analysis
    mag = 15; fps = 4;
    Nbins = 1;
    minLength = 75; % Min length of tracks
    Start = 1;      % Start Sequence
    DirEnd = 0;     % Value to subtract from length(total_directory)
    SwitchAxis = 0; % Switch x and y due to camera orientation (Switch = 1 .... Don't change = 0)
    micron_per_pix = 3.45; % FLIR

[BinnedTracks, ListYs, TracksLength] = TrackChannelDivideBin(total_directory, minLength,Nbins, mag, fps, Start, DirEnd, SwitchAxis,micron_per_pix);
'Done Binning'
% Calculate MSD
[MSD_x,MSD_x_adj,MSD_y,MSD_y_adj,MSD_XY,MSD_XY_adj,maxTime,CellSize,NumberDataPoints] = MSD_ViscosityCalc(BinnedTracks,TracksLength);
'Done Calc'

% Calculate Diffusion Coefficient
clc;
t = 1:maxTime; t_s = t./fps;
min = 1; max = 74;
kT = 4.11*10^-21;   % Thermal Energy [J=kg*m^2/s^2]
D = 0.5;            % diameter [microns]
a = 1.0*10^-6*D/2;  % radius [meters]
Diffusion = {};eta = {}; Diffusion_EP = {}; eta_EP = {};
for j = 1:length(MSD_y_adj)
    catMSDy = []; p = []; p_err = [];D = []; D_STD = [];
    catMSDy = 0.5*vertcat(MSD_y_adj{1,j});
    [p,p_err] = LinearFitError(t_s(min:max),(catMSDy(min:max))); % p(1) = Diffusion Coefficient [um^2/s]
    D = p(1) * (1.0*10^-6).^2; D_STD = p_err(1) * (1.0*10^-6).^2; % Diffusion & Prop. of Error [m^2/s]
    eta{1,j} = (kT./(6.*pi.*a.*D))/1E-3 % [cP]
    eta_EP{1,j} = sqrt((-kT./(6.*pi.*a) * (1/D)^2 * D_STD)^2) / 1.0E-3 % [cP]
    
    Diffusion{1,j} = D; % Diffusion [m^2/s]
    Diffusion_EP{1,j} = D_STD; % Propagation of Error: Diffusion [m^2/s]
end
% Plot the MSD vs time
close all;clc
min = 1; maxTime = 74;
t = 1:maxTime; t_s = t./fps;
PerPEOColor = [0 0 0; 
                0.635 0.078 0.184; 
                0.929 0.694 0.125; 
                0.466 0.674 0.188; 
                0.85 0.325 0.0980];
MarkerType = ['o','^','s','d','v'];
PlotMSD_Adjy = figure('color',[1 1 1]);
    p = zeros(1,CellSize(1))
    for j = 1:CellSize(1) % sequence to evaluate
        p(j) = plot(t_s(min:3:maxTime),(MSD_y_adj{1,j}(min:3:maxTime)),'LineStyle', 'none','marker',MarkerType(j),'color',PerPEOColor(j,:),...
            'LineWidth',1,'MarkerFaceColor',[1 1 1]);
        hold on;
    end

    fit = []; fit_y = [];  time_fit = [];
    for j = 1:CellSize(1) %sequence to evaluate
        fit(j,:) = polyfit(t_s(min:maxTime),MSD_y_adj{1,j}(min:maxTime),1);
        t_s_new = [t_s 22]
        y1 = polyval(fit(j,:),t_s_new);
        f = plot(t_s_new,y1,'-','color',PerPEOColor(j,:),'color',PerPEOColor(j,:),...
            'LineWidth',1);
        uistack(f,'bottom');
        time_fit = t_s_new; fit_y(:,j) = y1;
    end
    set(get(get(f(1:end),'Annotation'),'LegendInformation'),'IconDisplayStyle','off')
    lgd = legend(p(:),'0.0% PEO','0.1% PEO','0.25% PEO','0.5% PEO','0.82% PEO','location','northwest');
         
    hold off; box on;
    ax = gca; ax.FontSize = 7;
    xlim([-0.01 19]);ylim([0 33]);
    xticks([0 5 10 15]); yticks([0 5 10 15 20 25 30])
    set(gca, 'FontName', 'Arial')
    xlabel('Time (s)','FontSize',8);ylabel('MSD (\mum^2)','FontSize',8);
  
% Crop Figure
    alw = 1; % Border line width
    afs = 7; % Font Size

    lgd.EdgeColor = [1 1 1];
    lgd.FontName = 'Arial'; 
    lgd.FontSize = afs;
    
    markersize = 2;
    xlabel(''); ylabel('');
    mypapersize = [3.42 3]; % [width height]
    set(gca, 'FontName', 'Arial')
    set(gca,'LineWidth',alw);
    set(gca,'XColor',[0 0 0]);set(gca,'YColor',[0 0 0]);
    set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
    set(PlotMSD_Adjy,'Units','inches','Position',[0 0 mypapersize]);
% Save Figure
figdir = '.\ExtendedDataFigure4\'
print(PlotMSD_Adjy, '-painters','-dpdf', '-r1200', [figdir,'\MSDvsTime_V2.pdf']);
% Save Data
tab = 1;
dir_Excel = '.\Source data\ExtendedData\';
% Experiment
writematrix(t_s(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','A4')
writematrix(MSD_y_adj{1,1}(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','B4') % 0.0% PEO
writematrix(MSD_y_adj{1,2}(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','C4') % 0.1% PEO
writematrix(MSD_y_adj{1,3}(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','D4') % 0.25% PEO
writematrix(MSD_y_adj{1,4}(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','E4') % 0.5% PEO
writematrix(MSD_y_adj{1,5}(min:maxTime)',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','F4') % 0.82% PEO
% Fit
writematrix(time_fit',[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','H4')
writematrix(fit_y(:,1),[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','I4') % 0.0% PEO
writematrix(fit_y(:,2),[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','J4') % 0.1% PEO
writematrix(fit_y(:,3),[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','K4') % 0.25% PEO
writematrix(fit_y(:,4),[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','L4') % 0.5% PEO
writematrix(fit_y(:,5),[dir_Excel,'ExtDataFigure4.xlsx'],'Sheet',tab,'Range','M4') % 0.82% PEO
%% FUNCTION: 
function [BinnedTracks, ListYs, TracksLength] = TrackChannelDivideBin(total_directory, minLength,Nbins, mag, fps, Start, DirEnd, SwitchAxis,micron_per_pix)
% Divides the channel into Nbins number of bins and divides the tracks into
% each bin based on starting position
    count = 0;
    for i = Start:(length(total_directory)-DirEnd)
        count = count + 1;
        cd([total_directory(i).folder filesep total_directory(i).name]);
        load('ParticleTracks');
        for j = 1:length(tracks)
            if SwitchAxis == 1
            % Switching Axis because of Camera Orientation (Channel in 2500 pix direction)
                tracks(j).x_um = tracks(j).y.*(micron_per_pix/mag); tracks(j).y_um = tracks(j).x.*(micron_per_pix/mag);
                tracks(j).x_um_sub = tracks(j).x_um-tracks(j).x_um(1); tracks(j).y_um_sub = tracks(j).y_um-tracks(j).y_um(1);
                tracks(j).srootx2y2_um = sqrt((tracks(j).x_um.^2)+(tracks(j).y_um).^2); tracks(j).srootx2y2_um_sub = sqrt((tracks(j).x_um_sub.^2)+(tracks(j).y_um_sub).^2);
                TracksLength(j) = length(tracks(j).x_um);
            else
            % Camera orientation in ~2100 (smaller) direction
                tracks(j).x_um = tracks(j).x.*(micron_per_pix/mag); tracks(j).x_um_sub = tracks(j).x_um-tracks(j).x_um(1);
                tracks(j).y_um = tracks(j).y.*(micron_per_pix/mag); tracks(j).y_um_sub = tracks(j).y_um-tracks(j).y_um(1);
                tracks(j).srootx2y2_um = sqrt((tracks(j).x_um.^2)+(tracks(j).y_um).^2); tracks(j).srootx2y2_um_sub = sqrt((tracks(j).x_um_sub.^2)+(tracks(j).y_um_sub).^2);
                TracksLength(j) = length(tracks(j).x_um);
            end
        end
        
        Ys= cat(1,tracks(:).y_um); Ymin = min(Ys); Ymax = max(Ys);
        ListYs = linspace(Ymin,Ymax,Nbins+1); BinCount = 0;
        for k = 2:length(ListYs)
            BinCount = BinCount + 1;
            c = 0;
            for m = 1:length(tracks)
                if length(tracks(m).y_um) >= minLength && tracks(m).y_um(1) >= ListYs(k-1) && tracks(m).y_um(1) <= ListYs(k)
                    c = c + 1;
                    TracksBin(c).t = tracks(m).t;
                    TracksBin(c).x_um = tracks(m).x_um; TracksBin(c).x_um2 = (tracks(m).x_um).^2;
                    TracksBin(c).x_um_sub = tracks(m).x_um_sub; TracksBin(c).x_um_sub2 = (tracks(m).x_um_sub).^2;
                    TracksBin(c).y_um = tracks(m).y_um; TracksBin(c).y_um2 = (tracks(m).y_um).^2;
                    TracksBin(c).y_um_sub = tracks(m).y_um_sub; TracksBin(c).y_um_sub2 = (tracks(m).y_um_sub).^2;
                    TracksBin(c).magXY = tracks(m).srootx2y2_um; TracksBin(c).magXY2 = (tracks(m).srootx2y2_um).^2;
                    TracksBin(c).magXY_sub = tracks(m).srootx2y2_um_sub; TracksBin(c).magXY_sub2 = (tracks(m).srootx2y2_um_sub).^2;
                end
            end
            BinnedTracks{count,BinCount} = TracksBin;
            clear TracksBin
        end
    end
end
%% FUNCTION: Determine MSD
function [MSD_x,MSD_x_adj,MSD_y,MSD_y_adj,MSD_XY,MSD_XY_adj,maxTime,CellSize,NumberDataPoints] = MSD_ViscosityCalc(BinnedTracks,TracksLength)
% Uses circshift to take the difference between each cell by increments of
% delta_t, then it saves the values by increments of delta_t so they are
% stastically independent.
    maxTime = max(TracksLength);
    CellSize = size(BinnedTracks); % Dimension of the Cell; size(1,2) = [Row Column] = [Sequences Bins]
    x_sub_store = []; y_sub_store = []; magXY_sub_store = [];
    for i = 1:CellSize(1)     % Number of Sequences
        for g = 1:CellSize(2) % Number of Bins
            for t = 1:maxTime % Time intervals
                c = 0;
                for k = 1:length(BinnedTracks{i,g})
                            c = c + 1;
                            x_sub_store = BinnedTracks{i,g}(k).x_um -circshift(BinnedTracks{i,g}(k).x_um,t);
                                Delta_X{c} = x_sub_store(t+1:t:end);
                            
                             y_sub_store = BinnedTracks{i,g}(k).y_um -circshift(BinnedTracks{i,g}(k).y_um,t);
                                 Delta_Y{c} = y_sub_store(t+1:t:end);
                                 
                             magXY_sub_store = BinnedTracks{i,g}(k).magXY -circshift(BinnedTracks{i,g}(k).magXY,t);
                                 Delta_magXY{c} = magXY_sub_store(t+1:t:end);
                end
                                                
                x2_sub_bin(g,t) = mean(vertcat(Delta_X{:}).^2); x2_sub_bin_var(g,t) = var(vertcat(Delta_X{:}));
                y2_sub_bin(g,t) = mean(vertcat(Delta_Y{:}).^2); y2_sub_bin_var(g,t) =  var(vertcat(Delta_Y{:}));
                magXY2_sub_bin(g,t) = mean(vertcat(Delta_magXY{:}).^2); magXY2_bin_var(g,t) =  var(vertcat(Delta_magXY{:}));
                NumberDataPoint(g,t) = length(vertcat(Delta_X{:}));
                
                clear x_sub_store y_sub_store magXY_store Delta_X Delta_Y Delta_magXY
            end
        end
        MSD_x{i} = x2_sub_bin; MSD_x_adj{i} = x2_sub_bin_var;    
        MSD_y{i} = y2_sub_bin; MSD_y_adj{i} = y2_sub_bin_var;
        MSD_XY{i} = magXY2_sub_bin; MSD_XY_adj{i} = magXY2_bin_var;
        NumberDataPoints{i} = NumberDataPoint;
        clear x2_sub_bin x2_sub_bin_var y2_sub_bin y2_sub_bin_var magXY2_sub_bin magXY2_bin_var NumberDataPoint
    end
end

%% FUNCTION: 
function [p,p_err] = LinearFitError(x,y)
% This function calculates the error of the fit to a linear line: Y = BX + A
    N = length(x);
    delta = N*sum(x.^2)-sum(x)^2;

    A = (sum(x.^2)*sum(y) - sum(x)*sum(x.*y))/delta; % intercept 
    B = (N*sum(x.*y) - sum(x)*sum(y))/delta;         % slope
    p = [B A];                                       % mirror format of polyfit in matlab

    sigma_y = sqrt(1/(N-2) * sum((y-A-B*x).^2));

    sigma_A = sigma_y*sqrt(sum(x.^2)/delta); % intercept error
    sigma_B = sigma_y*sqrt(N/delta);         % slope error
    p_err = [sigma_B sigma_A];
end